Introduction à la modélisation statistique bayésienne

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

• Null Hypothesis Significance Testing
• Bayes Factor
• Comparaison de modèles : Overfitting
• Théorie de l'information
  • Incertitude
  • Divergence
  • Déviance
  • Régularisation
  • Critères d'information
• Application

On s'intéresse aux différences de taille entre hommes et femmes. On va mesurer 100 femmes et 100 hommes.

N = 100
men <- rnorm(N, 175, 10) # 100 tailles d'hommes
women <- rnorm(N, 170, 10) # 100 tailles de femmes

t.test(men, women)

    Welch Two Sample t-test

data:  men and women
t = 2.1475, df = 197.71, p-value = 0.03297
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2655634 6.2364837
sample estimates:
mean of x mean of y 
 173.0157  169.7646 

tibble(
    value = c(men, women),
    gender = c(rep("men", N), rep("women", N))
    ) %>%
    ggplot(aes(x = value, fill = gender) ) +
    geom_histogram() +
    theme_bw(base_size = 20)

On va simuler des t-valeurs issues de données générées selon l'hypothèse de non-différence entre hommes et femmes.

nSims <- 1e4 # number of simulations
t <- rep(NA, nSims) # initialising an empty vector

t <- replicate(
    nSims, 
    t.test(
        rnorm(N, 170, 10), 
        rnorm(N, 170, 10) 
        )$statistic
    )

La valeur d'un T-test est donnée par la formule :

\[ t = \frac{\mu_W - \mu_M}{S} \]

Où \( S \) est une fonction de la taille des échantillons et des variances des données.

data.frame(t = t) %>%
    ggplot(aes(x = t) ) +
    geom_histogram() +
    theme_bw(base_size = 25)

La distribution de Student t

Density, distribution function, quantile function and random generation for the t distribution with df degrees of freedom (and optional non-centrality parameter ncp).

Usage

dt(x, df, ncp, log = FALSE)
pt(q, df, ncp, lower.tail = TRUE, log.p = FALSE)
qt(p, df, ncp, lower.tail = TRUE, log.p = FALSE)
rt(n, df, ncp)

df_obs = t.test(men, women)$parameter
data.frame(t = t) %>%
    ggplot(aes(x = t)) +
    geom_histogram(aes(y = stat(density)), 
                   alpha = 0.4) +
    stat_function(fun = dt, 
        args = list(df = df_obs), size = 1.5) +
    theme_bw(base_size = 25)

tobs <- t.test(men, women)$statistic # observed t-value
tobs %>% as.numeric

[1] 2.147452

A p-value (in black) gives the probability of observing the data we observed (or more extreme data), given that the null hypothesis is true.

t.test(men, women)$p.value

[1] 0.03297456

tvalue <- abs(t.test(men, women)$statistic)
2 * integrate(dt, tvalue, Inf, df = df_obs)$value

[1] 0.03297456

2 * (1 - pt(abs(t.test(men, women)$statistic),df_obs) )

         t 
0.03297456 

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

• Null Hypothesis Significance Testing
• Bayes Factor
• Comparaison de modèles : Overfitting
• Théorie de l'information
  • Incertitude
  • Divergence
  • Déviance
  • Régularisation
  • Critères d'information
• Application

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

• Null Hypothesis Significance Testing
• Bayes Factor
• Comparaison de modèles : Overfitting
• Théorie de l'information
  • Incertitude
  • Divergence
  • Déviance
  • Régularisation
  • Critères d'information
• Application

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

• Null Hypothesis Significance Testing
• Bayes Factor
• Comparaison de modèles : Overfitting
• Théorie de l'information
  • Incertitude
  • Divergence
  • Déviance
  • Régularisation
  • Critères d'information
• Applica*tion

Courbe de divergence

t <- 
  tibble(
      p_1  = .3,
      p_2  = .7,
      q_1  = seq(from = .01, to = .99, by = .01)) %>% 
  mutate(q_2  = 1 - q_1) %>%
  mutate(d_kl = (p_1 * log(p_1 / q_1)) + (p_2 * log(p_2 / q_2)))

t %>% 
  ggplot(aes(x = q_1, y = d_kl)) +
  geom_vline(xintercept = .3, color = "red", linetype = 2) +
  geom_line(color = "steelblue", size = 1.5) +
  annotate(geom = "text", x = .4, y = 1.5, label = "q = p",
           color = "red", family = "Courier", size = 3.5) +
  labs(x = "q[1]",
       y = "Divergence of q from p")

Application

library(brms)
data(cars)

b <- 
  brm(data = cars, family = gaussian,
      dist ~ 1 + speed,
      prior = c(prior(normal(0, 100), class = Intercept),
                prior(normal(0, 10), class = b),
                prior(uniform(0, 30), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      seed = 6)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: dist ~ 1 + speed 
   Data: cars (Number of observations: 50) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept   -17.40      7.07   -30.97    -3.67       1920 1.00
speed         3.92      0.44     3.07     4.75       1783 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma    15.83      1.69    12.98    19.47       2068 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Calcul du nombre de paramètres effectifs

La deuxième partie du calcul du WAIC est le nombre de paramètres effectif, \( p_{WAIC} \). On définit \( V(y_{i}) \) comme la variance de la log-likelihood pour chaque observation \( i \) de l'échantillon d'entraînement.

\[ p_{\text{WAIC}} = \sum_{i=1}^{N} V(y_{i}) \]

dfvar <-
  ll %>%
  summarise_all(var) %>%
  gather(key, vars) %>%
  select(vars) 

pwaic <-
  dfvar %>%
  sum()

pwaic

[1] 3.418859

Calcul de WAIC

Ensuite, le WAIC est défini par:

\[ \text{WAIC} = -2(\text{lppd} - p_{\text{WAIC}}) \]

(WAIC <- -2 * (lppd - pwaic))

[1] 420.205

Le WAIC est aussi un estimateur de la déviance out-of-sample. La fonction waic du package brms permet de faire les calculs.

waic(b)

Computed from 4000 by 50 log-likelihood matrix

          Estimate   SE
elpd_waic   -210.1  6.4
p_waic         3.4  1.2
waic         420.2 12.7

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

• Null Hypothesis Significance Testing
• Bayes Factor
• Comparaison de modèles : Overfitting
• Théorie de l'information
  • Incertitude
  • Divergence
  • Déviance
  • Régularisation
  • Critères d'information
• Application

inits <- list(Intercept = mean(d$kcal.per.g),
              sigma     = sd(d$kcal.per.g))

inits_list <-list(inits, inits, inits, inits)

b6.11 <- 
  brm(data = d, family = gaussian,
      kcal.per.g ~ 1,
      prior = c(prior(uniform(-1000, 1000), class = Intercept),
                prior(uniform(0, 100), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      inits = inits_list,
      seed = 6)

inits <- list(Intercept = mean(d$kcal.per.g),
              neocortex = 0,
              sigma     = sd(d$kcal.per.g))
inits_list <-list(inits, inits, inits, inits)

b6.12 <- 
  brm(data = d, family = gaussian,
      kcal.per.g ~ 1 + neocortex,
      prior = c(prior(uniform(-1000, 1000), class = Intercept),
                prior(uniform(-1000, 1000), class = b),
                prior(uniform(0, 100), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      inits = inits_list,
      seed = 6)

inits <- list(Intercept   = mean(d$kcal.per.g),
              `log(mass)` = 0,
              sigma       = sd(d$kcal.per.g))
inits_list <-list(inits, inits, inits, inits)

b6.13 <-
  update(b6.12, 
         newdata = d,
         formula = kcal.per.g ~ 1 + log(mass),
         inits   = inits_list)

inits <- list(Intercept   = mean(d$kcal.per.g),
              neocortex   = 0,
              `log(mass)` = 0,
              sigma       = sd(d$kcal.per.g))
inits_list <-list(inits, inits, inits, inits)

b6.14 <- 
  update(b6.13, 
         newdata = d,
         formula = kcal.per.g ~ 1 + neocortex + log(mass),
         inits   = inits_list)

w[, 7:8] %>% 
  data.frame() %>% 
  rownames_to_column(var = "model_name") %>% 

  ggplot(aes(x    = model_name, 
             y    = waic, 
             ymin = waic - se_waic, 
             ymax = waic + se_waic)) +
  geom_pointrange(shape = 21, color = "steelblue") +
  coord_flip() +
  labs(x = NULL, y = NULL,
       title = "WAIC plot") +
  theme_classic(base_size  = 25) +
  theme(text             = element_text(family = "Courier"),
        axis.ticks.y     = element_blank(),
        panel.background = element_rect())

Le poids d'un modèle est une estimation de la probabilité que ce modèle fera les meilleures prédictions possibles sur un nouveau jeu de données, conditionnellement au set de modèles considéré.

\[ w_{i} = \frac{\exp (-\frac{1}{2}\text{dWAIC}_{i})}{\sum_{j=1}^{m} \exp (-\frac{1}{2}\text{dWAIC}_{j})} \]

model_weights(b6.11, b6.12, b6.13, b6.14, 
              weights = "waic") %>% 
  round(digits = 2)

b6.11 b6.12 b6.13 b6.14 
 0.02  0.01  0.02  0.96 

library(broom)
my_coef_tab <-
  tibble(model = c("b6.11", "b6.12", "b6.13", "b6.14")) %>% 
  mutate(fit   = purrr::map(model, get)) %>% 
  mutate(tidy  = purrr::map(fit, tidy)) %>% 
  unnest(tidy) %>% 
  filter(term != "lp__")

my_coef_tab %>%
  complete(term = distinct(., term), model) %>%
  select(model, term, estimate) %>%
  mutate(estimate = round(estimate, digits = 2)) %>%
  spread(key = model, value = estimate)

# A tibble: 4 x 5
  term        b6.11 b6.12  b6.13 b6.14
  <chr>       <dbl> <dbl>  <dbl> <dbl>
1 b_Intercept  0.66  0.35   0.71 -1.07
2 b_logmass   NA    NA     -0.03 -0.1 
3 b_neocortex NA     0.46  NA     2.77
4 sigma        0.19  0.19   0.18  0.14

Pourquoi ne conserver uniquement le premier modèle et oublier les autres ? Une autre stratégie consisterait à pondérer les prédictions des modèles par leurs poids respectifs. C'est ce qu'on appelle le moyennage de modèles (model averaging).

• Calculer le WAIC de chaque modèle
• Calculer le poids de chaque modèle
• Simuler des données à partir de chaque modèle
• Combiner ces valeurs simulées dans un ensemble de prédictions, en utilisant les poids des modèles comme proportions

# we need new data for both the `fitted()` and `pp_average()` functions
nd <- 
  tibble(neocortex = seq(from = .5, to = .8, length.out = 30),
         mass      = 4.5)

# we'll get the `b6.14`-implied trajectory with `fitted()`
f <-
  fitted(b6.14, newdata = nd) %>%
  as_tibble() %>%
  bind_cols(nd)

# the model-average trajectory comes from `pp_average()`
mdl_average <- pp_average(b6.11, b6.12, b6.13, b6.14,
           weights = "waic",
           method  = "fitted",  
           newdata = nd) %>%
    as_tibble() %>%
    bind_cols(nd) 

Voici les prédictions de tous les modèles considérés, pondérés par leur poids respectif. Comme le modèle b6.14 concentrait quasiment tout le poids, il fait sens que cette prédiction moyennée soit similaire aux prédictions du modèle b6.14.

mdl_average %>%
ggplot(aes(x = neocortex, y = Estimate)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), 
              fill  = "steelblue", alpha = 1/4) +
  geom_line(color   = "darkblue") +
  geom_ribbon(data  = f, aes(ymin = Q2.5, ymax = Q97.5),
              color = "darkgreen", fill = "green", alpha = .1, linetype = 2) +
  geom_line(data = f,
              color = "darkgreen", linetype = 2) +
  geom_point(data = d, aes(y = kcal.per.g), 
             size = 5, color = "orange") +
  labs(y = "kcal.per.g") +
  coord_cartesian(xlim = range(d$neocortex), 
                  ylim = range(d$kcal.per.g))

En bleu le modèle le modèle moyenné et en vert le modèle B6.14

Etude des donnés de Howell

1. Chargement des donnés
   Charger les donnés et créer deux jeux de donnés aléatoirement choisies de même taille.
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

Pour charger les donnés, utilisez la fonction

read.csv()

Pour selectionner les indice au hasard

sample()

Etude des données de Howell

1. Chargement des données
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

Howell1 <- read.csv("Data/Howell1.csv")

d <-
    Howell1 %>%
    mutate(age = scale(age) )

set.seed(1000)
i <- sample(1:nrow(d), size = nrow(d) / 2)

d1 <- d[i, ] # training sample
d2 <- d[-i, ] # test sample

Nous avons maintenant deux dataframes, de 272 lignes chacune. On va utiliser d1 pour fitter nos modèles et d2 pour les évaluer.

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
   Soit \( h_{i} \) les valeurs de taille et \( x_{i} \) les valeurs centrées d'âge, sur la ligne \( i \). Construisez les modèles suivants avec d1, en utilisant brm() et des priors faiblement régularisateurs.
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

\[ \begin{align} \mathcal{M_{1}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} \\ \mathcal{M_{2}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} \\ \mathcal{M_{3}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} \\ \mathcal{M_{4}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} x_{i}^{4} \\ \end{align} \]

Utiliser un intercept initial basé sur la moyenne des hauteurs, et une variance initiale basée sur la variance des hauteurs.

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
   Soit \( h_{i} \) les valeurs de taille et \( x_{i} \) les valeurs centrées d'âge, sur la ligne \( i \). Construisez les modèles suivants avec d1, en utilisant brm() et des priors faiblement régularisateurs.
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

inits <- list(Intercept = mean(d$height),
              age       = 0,
              sigma     = sd(d$height))

inits_list <-list(inits, inits, inits, inits)

mod_1 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age,
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list)

mod_2 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list
    )

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
   Soit \( h_{i} \) les valeurs de taille et \( x_{i} \) les valeurs centrées d'âge, sur la ligne \( i \). Construisez les modèles suivants avec d1, en utilisant brm() et des priors faiblement régularisateurs.
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

mod_3 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2) + I(age^3),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list)

mod_4 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2) + I(age^3) + I(age^4),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list
    )

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

Utiliser la méthode suivante pour calculer la valeur de WAIC et l'ajouter au modèle

model <- add_criterion(model, "waic")

Utiliser la méthod suivante pour effectuer la comparaison

loo_compare(models, criterion = "waic")

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

# compute and save the WAIC information for the next three models
mod_1 <- add_criterion(mod_1, "waic")
mod_2 <- add_criterion(mod_2, "waic")
mod_3 <- add_criterion(mod_3, "waic")
mod_4 <- add_criterion(mod_4, "waic")

# compare the WAIC estimates
w <- loo_compare(mod_1, mod_2, mod_3, mod_4,
                 criterion = "waic")

print(w)

      elpd_diff se_diff
mod_4    0.0       0.0 
mod_3  -40.3       9.3 
mod_2 -244.6      19.0 
mod_1 -509.0      21.0 

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

w[, 7:8] %>% 
  data.frame() %>% 
  rownames_to_column(var = "model_name") %>% 

  ggplot(aes(x = model_name, y = waic, 
             ymin = waic - se_waic, 
             ymax = waic + se_waic)) +
  geom_pointrange(shape = 21, color = "steelblue") +
  coord_flip() +
  labs(x = NULL, y = NULL, title = "WAIC plot") +
  theme_classic(base_size  = 25) +
  theme(text             = element_text(family = "Courier"),
        axis.ticks.y     = element_blank(),
        panel.background = element_rect())

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

Utiliser la méthode suivante pour calculer le poids respectif de chaque modèle

model_weights(models, weights = "waic")

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen

Le poids respectif de chaque modèle

model_weights(mod_1, mod_2, mod_3, mod_4, 
              weights = "waic") %>% 
  round(digits = 2)

mod_1 mod_2 mod_3 mod_4 
    0     0     0     1 

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen
   
   • Construire les modèles de degré 5 et 6
   • Calculer les poids respectifs de tous les modèles
   • Afficher le modèle moyen

mod_5 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2) + I(age^3) + I(age^4) + I(age^5),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list)

mod_6 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2) + I(age^3) + I(age^4) + I(age^5) + I(age^6),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list
    )

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen
   
   • Construire les modèles de degré 5 et 6
   • Calculer les poids respectifs de tous les modèles
   • Afficher le modèle moyen

# compute and save the WAIC information for the next three models
mod_5 <- add_criterion(mod_5, "waic")
mod_6 <- add_criterion(mod_6, "waic")

# compare the WAIC estimates
w <- loo_compare(mod_1, mod_2, mod_3, mod_4, mod_5, mod_6,
                 criterion = "waic")
print(w)

      elpd_diff se_diff
mod_6    0.0       0.0 
mod_4   -0.4       1.6 
mod_5   -0.9       1.5 
mod_3  -40.7       9.6 
mod_2 -245.0      18.9 
mod_1 -509.4      20.8 

model_weights(mod_1, mod_2, mod_3, mod_4, mod_5, mod_6, 
              weights = "waic") %>% 
  round(digits = 2)

mod_1 mod_2 mod_3 mod_4 mod_5 mod_6 
 0.00  0.00  0.00  0.32  0.20  0.48 

Etude des donnés de Howell

1. Chargement des donnés
2. Crée des modèles polynomiaux de degré 1 à 4
3. Comparer ces modèles en utilisant le WAIC.
4. Comparer les rangs des modèles et leurs poids.
5. Modèle moyen
   
   • Construire les modèles de degré 5 et 6
   • Calculer les poids respectifs de tous les modèles
   • Afficher le modèle moyen

# we need new data for both the `fitted()` and `pp_average()` functions
nd <- 
  tibble(age = seq(from = -2., to = 3., length.out = 50))

# we'll get the `b6.14`-implied trajectory with `fitted()`
f <-
  fitted(mod_6, newdata = nd) %>%
  as_tibble() %>%
  bind_cols(nd)

# the model-average trajectory comes from `pp_average()`
mdl_average <- pp_average(mod_1, mod_2, mod_3, mod_4, mod_5, mod_6,
           weights = "waic",
           method  = "fitted",  
           newdata = nd) %>%
    as_tibble() %>%
    bind_cols(nd) 

Voici les prédictions de tous les modèles considérés, pondérés par leur poids respectif. Comme le modèle b6.14 concentrait quasiment tout le poids, il fait sens que cette prédiction moyennée soit similaire aux prédictions du modèle b6.14.

mdl_average %>%
ggplot(aes(x = age, y = Estimate)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), 
              fill  = "steelblue", alpha = 1/4) +
  geom_line(color   = "darkblue") +
  geom_ribbon(data  = f, aes(ymin = Q2.5, ymax = Q97.5),
              color = "darkgreen", fill = "green", alpha = .1, linetype = 2) +
  geom_line(data = f,
              color = "darkgreen", linetype = 2) +
  geom_point(data = d, aes(y = height), 
             size = 2, color = "orange") +
  labs(y = "kcal.per.g") +
  coord_cartesian(xlim = range(d$age), 
                  ylim = range(d$height)) +